Modeling Permafrost: Soil, Ice, and Some Really Hard Mathematics

By Naren Vohra and Malgorzata Peszynska

Permafrost—i.e., ground that remains frozen for two or more years—is abundant in the Arctic and covers 85 percent of Alaska. It is a highly heterogeneous, multicomponent system that consists of soil layers (silt, clay, peat, and so forth) and ground ice in the form of ice lenses or massive ice wedges on the order of one meter (m). Permafrost thaws when ground temperatures increase, leading to ground subsidence and deformation that can damage infrastructure—such as buildings, railways, and pipelines (see Figure 1a)—and erode natural landscape features like pingos, which become large marshy lakes called thermokarsts.

Efficient computational models for permafrost thaw require accurate mathematical representations of the different physical processes and their intricate coupling across many spatial and temporal scales. These processes include heat conduction with phase change (Tp), hydrological flow (H), and mechanical deformation (M), which all have a strong interdependence that leads to the fully coupled thermo-hydro-mechanical (TpHM) system (see Figure 1b). While TpHM solvers do exist in the commercial packages that the geotechnical community utilizes for case-specific engineering projects, they are not easily accessible to those who are interested in permafrost-specific simulation scenarios across a range of scales. Additionally, the current thermo-hydrological (TpH) models do not address questions of robustness, accuracy, and algorithmic stability.

Using a wealth of applied and computational mathematical techniques, we aimed to develop and analyze numerical schemes with provable convergence and robustness properties that can handle a range of simulation scenarios. These schemes then enable us to quantify the tradeoff between accuracy and efficiency depending on discretization and the spatial and temporal scales, thereby allowing faster large-scale simulations and, eventually, the massive use of experimental data. 

Figure 1. Permafrost thaw disturbs local infrastructure and natural landscape features. 1a. Permafrost thaw due to improper road construction leading to the development of a thermokarst. 1b. The different processes and relationships in the thermo-hydro-mechanical coupling that governs permafrost thaw. Figure 1a courtesy of John Cloud through the National Oceanic and Atmospheric Administration Photo Library, and Figure 1b courtesy of the authors.

Details of the Model 

Beginning with Tp, we can model heat conduction with advection and phase transitions via the nonlinear degenerate parabolic system

\[\partial_tw-\nabla \cdot (k\nabla\theta)+\nabla \cdot (c\theta q_f)=0, \enspace w \in \alpha(\theta).\tag1\]

Here, \(\theta\) is the temperature (in \(^\circ\) Celsius (\(\textrm{C}\))), \(w\) is the enthalpy per unit volume (in \(\textrm{J/m}^3\)), \(k\) is the thermal conductivity (in \(\textrm{J/m s } ^\circ \textrm{C}\)), \(c\) is the volumetric heat capacity of liquid water (in \(\textrm{J/m}^3\) \(^\circ \textrm{C}\)), \(q_f\) is the hydrological flux (in \(\textrm{m/s}\)), and \(\alpha\) is the relationship between temperature and enthalpy. When \(q_f=0\) in \((1)\), the system reduces to the Tp process [1, 3], which is mathematically well-studied in bulk water [8].

To obtain \(q_f\), we model hydrological flow and mechanical deformation (HM) with the well-known Biot’s poroelasticity equations, which comprise a strongly coupled system of linear elasticity and compressible Darcy flow in saturated porous media:

\[-\nabla \cdot [\lambda (\nabla \cdot u)I+2\mu\epsilon(u)]+\alpha_B\nabla p=0,\tag2\]

\[\partial_t(c_0p+\alpha_B\nabla \cdot u)+\nabla \cdot q_f=0.\tag3\]

Here, \(u\) is the displacement (in \(\textrm{m}\)), \(p\) is the pressure (in pascals (\(\textrm{Pa}\))), \(\lambda\) and \(\mu\) are the Lamé coefficients, \(\epsilon(u)\) is the linearized strain tensor, \(\alpha_B\) is the Biot-Willis constant, \(c_0\) is the specific storage coefficient (in \(1/\textrm{Pa}\)), and \(I\) is the identity matrix. The flux is given by \(q_f=-\kappa\mu_f^{-1}(\nabla p)\) via Darcy’s law, where \(\kappa\) is the permeability (in \(\textrm{m}^2\)) of the media and \(\mu_f\) is the viscosity (in \(\textrm{Pa s}\)) of liquid water.

The fully coupled TpHM system for permafrost introduces suitable constitutive relationships through additional coupling terms in \((1)\), \((2)\), and \((3)\). These relationships include hydrological phenomena such as cryosuction (the movement of water towards an ice front), interdependence of hydromechanical parameters (like permeability and Lamé coefficients) on temperature and water content, and the deterioration of these parameters during ground ice thaw. Such phenomena are evident from lab experiments and field data or identifiable during the upscaling of processes at pore-scale, starting with X-ray computed tomography images [5].

Once we define appropriate constitutive relationships, we can then find a numerical approximation. Previous studies have utilized piecewise-linear P1 Galerkin finite elements [3, 4], such as for the Stefan problem in bulk water [7]. However, heterogeneity and the multiphysics nature of TpHM processes in permafrost pose several challenges, which are best served with conservative approaches that are robust in incompressible limits or when \(\kappa \rightarrow 0\) [9]. These complications add to the challenges that arise directly from the constitutive relationships in Tp — such as the nonlinear relationship \(\alpha(\theta)\) in \((1)\), which is non-smooth and multivalued at material interfaces. The typical approach to regularizing \(\alpha\) introduces inconsistency and modeling errors that propagate in a manner that is difficult to track in a coupled system.

Figure 2. The effect of a warming climate on permafrost thaw. We consider a heterogeneous domain that includes an ice wedge, and calculate the temperature and enthalpy profiles at the time of maximum thaw depth (indicated by the dashed lines in the enthalpy plot). Figure courtesy of the authors.

Model Results

To address the aforementioned challenges, we focused on mixed finite elements for TpHM. Beginning with Tp, we developed a robust mixed finite element scheme that seeks temperature and enthalpy in the space of piecewise constants \(\textrm{P}0\) [1], with fluxes in the lowest-order Raviart-Thomas space \(\textrm{RT}_{[0]}\). Our \(\textrm{P}0\)-\(\textrm{P}0\) scheme is conservative, can accurately handle phase changes and heterogeneity, and can function as a cell-centered finite difference scheme. In Figure 2, we simulate Tp with realistic soil parameters to find the thickness of the non-isothermal portion of frozen clay (called the active layer) under an extreme climate warming scenario. Our results predict that the maximum thaw depth increases by 40 percent after two years if the temperature increases at a rate of \(5 \ ^\circ \textrm{C/year}\).

Moving on to HM, we discretize \((2)\) with a three-field formulation that utilizes \(\textrm{Q}1\)-\(\textrm{P}0\)-\(\textrm{RT}_{[0]}\) elements [6], where \(\textrm{Q}1\) is the space of continuous bilinear polynomials. Figure 3 illustrates a thawing scenario for a \(20 \times 20 \times 40 \textrm{ m}^3\) column of frozen silt with a (somewhat high) thaw rate of \(24 \textrm{ m/year}\) and a large external stress of \(1\) megapascal (a very high value, say, for a small house). Using these values, our model predicts subsidence of approximately \(4 \textrm{ m}\) at the end of one year.

Figure 3. Subsidence due to permafrost thaw under external stress, simulated with our poroelasticity code (based on deal.ii) and textured with Blender for interpretation. Figure courtesy of the authors.

This demonstration serves two purposes. The first is positive: it shows robustness of the three-field formulation in the context of heterogeneity and relevant physical parameters. The second, however, is somewhat negative: the demonstration reveals the inability of linear poroelasticity to reproduce subsidence of the level that is present in permafrost regions (on the order of 10 centimeters) unless we apply an unrealistically high magnitude of input that is typically absent in these areas. We hypothesize that one must consider additional model elements—such as full TpHM coupling, cryosuction, water flow, and nonlinear extensions of \((2)\)—and account for physical heterogeneity to avoid this issue.

In particular, we recognize the need to better understand and model massive ice wedges and ice lenses that become flowing water after thawing and destabilize the soil’s mechanical response. However, this factor leads to further questions. We can easily simulate the Tp process for ice wedges (see Figure 4 for an illustration with our \(\textrm{P}0\)-\(\textrm{P}0\) scheme), but what happens after the ice melts? Does the liquid water pond on the surface? How can we model the associated subsidence? In short, how should we incorporate the ice wedges and their fate into the HM framework? And how can we accurately and robustly simulate the coupled processes [2]? 

Figure 4. An ice wedge (outlined in black) thaws over the course of 10 years due to a ground surface temperature of \(20 \ ^\circ \textrm{C}\). A major challenge for thermo-hydro-mechanical solvers is modeling the excess liquid water and the subsidence due to thawing ground ice. Figure courtesy of the authors.

Conclusions and Future Work

We are currently working on algorithms that robustly couple our \(\textrm{P}0\)-\(\textrm{P}0\) Tp scheme with the \(\textrm{Q}1\)-\(\textrm{P}0\)-\(\textrm{RT}_{[0]}\) HM scheme to simulate TpHM coupling. We have demonstrated our \(\textrm{P}0\)-\(\textrm{P}0\) Tp algorithm’s ability to robustly handle ice wedges—even when coupled with flow TpH—and are now investigating its behavior in the HM framework.

In addition, we are collaborating with Elchin Jafarov and Brendan Rogers of Woodwell Climate Research Center to create a surrogate model of large-scale thaw settlement that is correlated to active layer depth. We will be able to use data from the Arctic to explore a data-to-data model of ground subsidence.

To conclude, permafrost modeling is a complex enterprise that involves well-understood individual components: the Tp, H, M, and HM processes. When confronted with data and reality, however, TpHM coupling presents unforeseen challenges to computational models and solvers. Nevertheless, these challenges inspire respect and awe within the geophysical and geotechnical communities, which continue to engage at the forefront of Arctic permafrost research.

SIAM Student Travel Awards are made possible in part by the generous support of our community. To make a gift to the Student Travel Fund, visit the SIAM website. 

Acknowledgments: This research was partially supported by NSF DMS-1912938: “Modeling with Constraints and Phase Transitions in Porous Media” and NSF DMS-1522734: “Phase transitions in porous media across multiple scales” under principal investigator Malgorzata Peszynska of Oregon State University. The research was also partially supported by funding catalyzed through the Audacious Project. We would also like to acknowledge Lisa Bigler of Sandia National Laboratories for her contributions to the early part of this research [1].

References
[1] Bigler, L., Peszynska, M., & Vohra, N. (2022). Heterogeneous Stefan problem and permafrost models with P0-P0 finite elements and fully implicit monolithic solver. Electron. Res. Arch., 30(4), 1477-1531.
[2] Brun, M.K., Ahmed, E., Berre, I., Nordbotten, J.M., & Radu, F.A. (2020). Monolithic and splitting solution schemes for fully coupled quasi-static thermo-poroelasticity with nonlinear convective transport. Comput. Math. Appl., 80(8), 1964-1984.
[3] Nicolsky, D.J., Romanovsky, V.E., & Tipenko, G.S. (2007). Using in-situ temperature measurements to estimate saturated soil thermal properties by solving a sequence of optimization problems. Cryosphere, 1(1), 41-58.
[4] Nicolsky, D.J., Romanovsky, V.E., Tipenko, G.S., & Walker, D.A. (2008). Modeling biogeophysical interactions in nonsorted circles in the Low Arctic. J. Geophys. Res. Biogeosci., 113(G3), G03S05.
[5] Peszynska, M., Trykozko, A., Iltis, G., Schlueter, S., & Wildenschild, D. (2016). Biofilm growth in porous media: Experiments, computational modeling at the porescale, and upscaling. Adv. Water Resour., 95(2), 288-301.
[6] Phillips, P.J., & Wheeler, M.F. (2007). A coupling of mixed and continuous Galerkin finite element methods for poroelasticity I: The continuous in time case. Comput. Geosci., 11(2), 131-144.
[7] Verdi, C., & Visintin, A. (1988). Error estimates for a semi-explicit numerical scheme for Stefan-type problems. Numer. Math., 52(2), 165-185.
[8] Visintin, A. (1996). Models of phase transitions (1st ed.). In Progress in nonlinear differential equations and their applications (Vol. 28). Boston, MA: Birkhäuser.
[9] Yi, S.-Y. (2017). A study of two modes of locking in poroelasticity. SIAM J. Numer. Anal., 55(4), 1915-1936.

Naren Vohra is a Ph.D. candidate in the Department of Mathematics at Oregon State University (OSU). His research interests include mathematical and computational modeling of multiphysics multiscale processes in porous media with applications to permafrost thaw. More information about his work is available on his website.

Malgorzata Peszynska is a professor of mathematics at OSU. She received her Ph.D. in 1992 from the University of Augsburg in Germany and has held positions at the Polish Academy of Sciences, Warsaw University of Technology, Purdue University, University of Texas at Austin, and National Science Foundation. She is a 2020 American Association for the Advancement of Science Honorary Fellow, the recipient of the 2021 SIAM Activity Group on Geosciences Prize and 2022 OSU Champion of Science Dean’s Award, and the 2022-2024 Joel Davis Faculty Scholar at OSU. She was also a SIAM honoree for Women's History Month in 2023. More information about her research, publications, and students can be found on her website.